function diff_im = anisodiff_unit(im, num_iter, delta_t, kappa, option)
% diff_im = anisodiff_unit(im, num_iter, delta_t, kappa, option)
%   conventional anisotropic diffusion (Perona & Malik) on a gray scale
%   image. A 2D network structure of 8 neighboring nodes is considered for 
%   diffusion conduction.
% 
%       ARGUMENT DESCRIPTION:
%              IM       - gray scale image (MxN).
%              NUM_ITER - number of iterations. 
%              DELTA_T  - integration constant (0 <= delta_t <= 1/7).
%                         Usually, due to numerical stability this 
%                         parameter is set to its maximum value.
%              KAPPA    - gradient modulus threshold that controls the conduction.
%              OPTION   - conduction coefficient functions proposed by Perona & Malik:
%                         1 - c(x,y,t) = exp(-(nablaI/kappa).^2),
%                             privileges high-contrast edges over low-contrast ones. 
%                         2 - c(x,y,t) = 1./(1 + (nablaI/kappa).^2),
%                             privileges wide regions over smaller ones. 
% 
%       OUTPUT DESCRIPTION:
%               DIFF_IM - (diffused) image with the largest scale-space parameter.
% 
%   Example
%   -------------
%   s = phantom(512) + randn(512);
%   num_iter = 15;
%   delta_t = 1/7;
%   kappa = 30;
%   option = 2;
%   ad = anisodiff2D(s,num_iter,delta_t,kappa,option);
%   figure, subplot 121, imshow(s,[]), subplot 122, imshow(ad,[])
%
% References: 
%   P. Perona and J. Malik. 
%   Scale-Space and Edge Detection Using Anisotropic Diffusion.
%   IEEE Transactions on Pattern Analysis and Machine Intelligence, 
%   12(7):629-639, July 1990.
% 
%   G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
%   Nonlinear Anisotropic Filtering of MRI Data.
%   IEEE Transactions on Medical Imaging,
%   11(2):221-232, June 1992.
%  
%   Codes adapted from Daniel Simoes Lopes
%   ICIST
%   Instituto Superior Tecnico - Universidade Tecnica de Lisboa
%   danlopes (at) civil ist utl pt
%   http://www.civil.ist.utl.pt/~danlopes
% 
%   Jinghao Lu 01/06/2016

% Convert input image to double.
im = double(im);

% PDE (partial differential equation) initial condition.
diff_im = im;

% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);

% 2D convolution masks - finite differences.
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];

% Anisotropic diffusion.
for t = 1:num_iter

        % Finite differences. [imfilter(.,.,'conv') can be replaced by conv2(.,.,'same')]
%         nablaN = imfilter(diff_im,hN,'conv');
%         nablaS = imfilter(diff_im,hS,'conv');   
%         nablaW = imfilter(diff_im,hW,'conv');
%         nablaE = imfilter(diff_im,hE,'conv');   
%         nablaNE = imfilter(diff_im,hNE,'conv');
%         nablaSE = imfilter(diff_im,hSE,'conv');   
%         nablaSW = imfilter(diff_im,hSW,'conv');
%         nablaNW = imfilter(diff_im,hNW,'conv'); 
        nablaN = conv2(diff_im, hN, 'same');
        nablaS = conv2(diff_im, hS, 'same');   
        nablaW = conv2(diff_im, hW, 'same');
        nablaE = conv2(diff_im, hE, 'same');   
        nablaNE = conv2(diff_im, hNE, 'same');
        nablaSE = conv2(diff_im, hSE, 'same');   
        nablaSW = conv2(diff_im, hSW,'same');
        nablaNW = conv2(diff_im, hNW, 'same'); 
        
        % Diffusion function.
        if option == 1
            cN = exp(-(nablaN/kappa).^2);
            cS = exp(-(nablaS/kappa).^2);
            cW = exp(-(nablaW/kappa).^2);
            cE = exp(-(nablaE/kappa).^2);
            cNE = exp(-(nablaNE/kappa).^2);
            cSE = exp(-(nablaSE/kappa).^2);
            cSW = exp(-(nablaSW/kappa).^2);
            cNW = exp(-(nablaNW/kappa).^2);
        elseif option == 2
            cN = 1./(1 + (nablaN/kappa).^2);
            cS = 1./(1 + (nablaS/kappa).^2);
            cW = 1./(1 + (nablaW/kappa).^2);
            cE = 1./(1 + (nablaE/kappa).^2);
            cNE = 1./(1 + (nablaNE/kappa).^2);
            cSE = 1./(1 + (nablaSE/kappa).^2);
            cSW = 1./(1 + (nablaSW/kappa).^2);
            cNW = 1./(1 + (nablaNW/kappa).^2);
        end

        % Discrete PDE solution.
        diff_im = diff_im + ...
                  delta_t*(...
                  (1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
                  (1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...
                  (1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...
                  (1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );
%                   delta_t*(...
%                   (1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
%                   (1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE);
           
        % Iteration warning.
%         fprintf('\rIteration %d\n',t);
end